package hr.irb.goTermList;

import hr.irb.geneOntology.*;
import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collection;
import java.util.Collections;
import java.util.Formatter;
import java.util.HashMap;
import java.util.HashSet;
import java.util.Iterator;
import java.util.LinkedHashSet;
import java.util.List;
import java.util.Locale;
import java.util.Map;
import java.util.Set;
import weka.filters.unsupervised.attribute.PrincipalComponents;
import weka.core.Attribute;
import weka.core.FastVector;
import weka.core.Instance;
import weka.core.Instances;
import weka.core.WekaException;
import weka.filters.Filter;
import weka.core.converters.CSVSaver;


/**
 * Takes a list of GO terms and prepares them for visualization in 1-, 2- or 3-
 * dimensional space, where coordinates are derived from term semantic
 * distances summarized by principal component analysis (PCA).
 * 
 * First, the pairwise distances of all listed GO terms are determined -
 * currently, the XML-RPC interface of FunSimMat service is used 
 * ( http://funsimmat.bioinf.mpi-inf.mpg.de/ ). In the future, the distances
 * might also be read from a local copy of the FSST database, or computed 
 * using GO4J.
 * 
 * Then the myInst is stored into an Instances object (from Weka) and principal
 * component analysis (PCA) is performed, retaining the first three PCs.
 * Alternatively, other dimensionality reduction methods (e.g. SOM, MDS)
 * Whatever myInst is stored as "properties" in each GO term is added to the
 * Instances afterwards (can be used in visualization - disc size, color etc.)
 * 
 * Finally, some GO terms are indicated as 'disposable'. This is determined by 
 * recursively eliminating one term from a pair of currently closest GO terms 
 * in the set, until the distance of the closest pair drops below a given 
 * threshold. Which GO term is eliminated from a pair is determined by checking
 * the value of a term's property, if exists (if not, one is removed at random).
 * (Future plan: if no property is supplied, the term with the smaller minimum
 * distance to all other terms in list is kept.)
 * 
 * @author Fran Supek (fran.supek[AT]irb.hr)
 */

public class TermListVisualizer {

  /** Some string constants defining properties of test results. */
  public final static String enrichString = "enrichment";
  public final static String log10SizeString = "log_size";
  public final static String pValueString = "log_pValue";
  public final static String flagFreqString = "flagFreq";
  public final static String inGroupFlaggedCntString = "inGroupFlagged";
  public final static String inGroupNoFlagCntString = "inGroupNoFlag";
  public final static String outOfGroupFlaggedCntString = "outOfGroupFlagged";
  public final static String outOfGroupNoFlagCntString = "outOfGroupNoFlag";
  public final static String isTestBString = "isTestB";

  public GoTermList termList;
  
  protected GeneOntology myGo;
  protected String description;
  protected String orgsInTotal;
  protected GoNamespace namespace;

  /** 
   * NOT SUPPORTED ANYMORE - can be used to compute GO Term semantic
   * similarities. (the code is inside & commented, it may be brought to life
   * on demand
   */
  // protected final GenePlex genePlex;
  
  /** Needed to compute the semantic similarities. */
  protected final GoTermSizes termSizes;
     
  /**
   * Requires that a List of GOTerms is provided, as
   * well as the GeneOntology instance they were taken from. <p>
   *
   * The TermListVisualizer will use the provided GoTermSizesObject to
   * calculate semantic distances. <p>
   * 
   * @param termList A list of selected GO terms.
   * @param go A GeneOntology object with all GO terms.
   * @param termSizes A GoTermSizes object necessary to compute semantic
   * distances.
   * @param desc A free-text description of the visualization. Used for 
   * header in output arff files, and in reporting of errors. Optional.
   * @param orgsInTotal I don't remember what this was for; it isn't used for
   * any calculations, it is just output as-is to Weka Instances.
   */
  public TermListVisualizer (GoTermList termList,
          GeneOntology go, GoTermSizes termSizes,
          String desc, String orgsInTotal) {
    this.termList = termList;
    myGo = go;
    if (desc==null)
      description = "";
    else
      description = desc;
    namespace = GoNamespace.findNamespaceOfList(termList);
    this.orgsInTotal = orgsInTotal;

    this.termSizes = termSizes;
  }  

  public TermListVisualizer (GoTermList termList,
          GeneOntology go, GoTermSizes termSizes) {
    this(termList, go, termSizes, "", "");
  }  



  protected Instances scatterPoints(Instances parInst, int xColumn, int yColumn)
          throws Exception {
      Instances myInst = parInst;
      int numInstances = myInst.numInstances();

      double[][] vectors = new double[numInstances][2];
      double[][] oldVectors = new double[numInstances][2];
      double dt2 = 0.1,                 // time differential used in "simulation"
              distanceThreshold = 0.04, // vectors distanced above this are not
              scaling = 0.8,            // scaling of the exponential damping (potential function)
              expFactor = 20;           // correction of the exponential factor in damping

      for (int i=0;i<numInstances;++i){
          vectors[i][0] = myInst.instance(i).value(xColumn);
          vectors[i][1] = myInst.instance(i).value(yColumn);
      }

      for (int counter = 0; counter <17; counter++){
          //newVectors = vectors, deep copy
          for (int c=0;c<vectors.length;++c){
              oldVectors[c][0] = vectors[c][0];
              oldVectors[c][1] = vectors[c][1];
          }
          for (int i=0;i<numInstances-1;++i){
              double[] point1 = oldVectors[i];
              for (int j=i+1;j<numInstances;++j){
                  double[] point2 = oldVectors[j];
                  double[] vectorDiff = {point2[0]-point1[0],point2[1]-point1[1]};
                  // distance between vectors
                  double distance = Math.sqrt(vectorDiff[0]*vectorDiff[0] + vectorDiff[1]*vectorDiff[1]);
                  if (distance==0)
                      distance=0.00001;
                  if (distance>distanceThreshold)
                      continue;         // skipping points beyond a certain distance
                  double accelerationValue = scaling * Math.exp(-expFactor*distance);
                  // acceleration vector
                  double[] acceleration = {accelerationValue*vectorDiff[0]/distance,accelerationValue*vectorDiff[1]/distance};
                  // move current vector by a/2 * dt2
                  vectors[i][0] -= 0.5*dt2*acceleration[0];
                  vectors[i][1] -= 0.5*dt2*acceleration[1];
                  // also, move the other vector too
                  vectors[j][0] += 0.5*dt2*acceleration[0];
                  vectors[j][1] += 0.5*dt2*acceleration[1];
              }
          }
      }

      for (int i=0;i<numInstances;++i){
          myInst.instance(i).setValue(xColumn, vectors[i][0]);
          myInst.instance(i).setValue(yColumn, vectors[i][1]);
      }
      return myInst;
  }
  
  
  
  /**
   * Reduces dimensionality of the supplied Instances object, currently by 
   * using Weka's Principal Component Analysis and by retaining at most first
   * numAxes PCs, regardless of the total % of information retained.
   */
  protected Instances reduceDimensionality(Instances parInst, int numAxes)
          throws Exception {

    /* Standardize myStand = new Standardize();
    myStand.setInputFormat(parInst);
    Instances myInst = Filter.useFilter(parInst, myStand); */
    
    //System.err.println("Started reduceDimensionality().");

    Instances myInst = parInst;
    myInst.setClassIndex(-1);

    PrincipalComponents myPCA = new PrincipalComponents();
    myPCA.setMaximumAttributes(numAxes); myPCA.setNormalize(false);
    myPCA.setMaximumAttributeNames(numAxes); myPCA.setInputFormat(myInst);
    myPCA.setVarianceCovered(0.999);    
    myInst = Filter.useFilter(myInst, myPCA);

    //System.err.println("Finished reduceDimensionality().");

    return myInst;
    
  }

  /**
   * Reduces dimensionality of the supplied Instances object by using MDSJ.
   */
  protected Instances reduceDimensionalityMds(Instances parInst, int numAxes)
          throws Exception {

    /* Standardize myStand = new Standardize();
    myStand.setInputFormat(parInst);
    Instances myInst = Filter.useFilter(parInst, myStand); */

    //System.err.println("Started reduceDimensionality().");

    double[][] data = new double[parInst.numInstances()][];
    for (int i = 0; i< parInst.numInstances(); i++) {
      data[i] = parInst.instance(i).toDoubleArray();
      for ( int j = 0; j<data[i].length; j++ ) {
        if (i==j) {
          data[i][j] = 0.0;
        } else {
          //data[i][j] = 1.0 - data[i][j];
          data[i][j] = 1.0 / Math.max( data[i][j], 0.1 );
        }
      }
    }

    double[][] newData = mdsj.MDSJ.stressMinimization(data, numAxes);

    FastVector atts = new FastVector();
    for (int j = 0; j < numAxes; j++) {
      Attribute curAtt = new Attribute(
              String.format( "axis_%02d", j+1 ) );
      atts.addElement(curAtt);
    }
    Instances myInst = new Instances("MDS_scaled", atts, parInst.numInstances());
    myInst.setClassIndex(-1);

    double[] vals;
    for (int i = 0; i< parInst.numInstances(); i++) {
      vals = new double[numAxes];
      for (int j = 0; j < numAxes; j++) {
        vals[j] = newData[j][i];
      }
      myInst.add( new Instance(1.0, vals) );
    }

    return myInst;

  }


  /**
   * Two values (x and y) of the supplied property are considered equal
   * if abs(x-y)/((x+y)/2) is less than the confidenceInterval.
   * It's recommended to use "absLogProp" on p-values if you want this to work
   * correctly.
   */
  final double confidenceInterval = 0.10;  // 10 % of average value

  /**
   * A genegroup is considered too general if it contains more than this
   * proportion of total genes. Such genegroups are considered to be
   * uninteresting.
   */
  final double genegroupGeneralityThreshold = 0.05;  // suggestions: 1.5 % or 5 %



  /**
   * Determines 'dispensability' of each GO term in the list. This is done by
   * recursively eliminating one term from a pair of currently most similar GO
   * terms in the set, and remembering the similarity of the pair as the
   * 'dispensability' of the eliminated term.
   *
   * Which GO term is eliminated from a pair is determined by
   * checking the value of a term's property, if exists (if not, one is removed
   * at random).
   *
   * @param termProperties - Pre-populate this with the property of interest.
   * @param criterion - Name of the GOTerm's property to check. Larger
   * values are assumed to be more desirable.
   * @param altCriterion - If two terms are similar by the criterion above,
   * check value of the altCriterion. You may pass null here.
   * @param parSimilMat - The distance matrix of GOTerms.
   * @param goTermSizes - Helps in deciding which terms to keep by avoiding
   * very general terms, and it also influences the outcome if terms are linked
   * by parent-child relationships. You may pass "null" here.*
   */
  protected void markDispensableTerms( GoTermProperties termProperties,
          String criterion, String altCriterion,
          TermSimilarityMatrix parSimilMat,
          GoTermSizes goTermSizes ){
    markDispensableTerms(termProperties, criterion, altCriterion, true,
            false, parSimilMat, goTermSizes);
  }


  /**
   * Determines 'dispensability' of each GO term in the list. This is done by 
   * recursively eliminating one term from a pair of currently most similar GO
   * terms in the set, and remembering the similarity of the pair as the
   * 'dispensability' of the eliminated term.
   * 
   * Which GO term is eliminated from a pair is determined by
   * checking the value of a term's property, if exists (if not, one is removed
   * at random).
   * 
   * @param criterion - Name of the GOTerm's property to check.
   * @param altCriterion - If two terms are similar by the criterion above,
   * check value of the altCriterion. You may pass null here.
   * @param keepGreater - If true, the GOTerm with a larger value of the
   * property wins (the property is 'desirable'). 
   * @param absLogProp - If true, transforms property into absolute value of
   * log(property) before comparing.
   * @param parSimilMat - The distance matrix of GOTerms.
   * @param goTermSizes - Helps in deciding which terms to keep by avoiding
   * very general terms, and it also influences the outcome if terms are linked
   * by parent-child relationships. You may pass "null" here.
   */
  protected void markDispensableTerms( GoTermProperties termProperties,
          String criterion, String altCriterion, Boolean keepGreater,
          Boolean absLogProp, TermSimilarityMatrix parSimilMat,
          GoTermSizes goTermSizes ){

    // this is to keep the original matrix from getting destroyed
    TermSimilarityMatrix similMat = parSimilMat.clone();

    double maxRemainingSimil = 0;
    int maxX, maxY;
            
    // set "dispensability" of all terms to 0, because otherwise one term
    // could be left without a dispensability value
    for ( GOTerm curTerm : termList ) {
      termProperties.put(curTerm, "dispensability", 0.0);
      termProperties.put(curTerm, "dispensedBy", 0.0);
    }


    /* If two GOTerms have opposite signs (ie. one is + an the other -) of
     * the property, set their pairwise correlation to 0.
     * Normally they shouldn't have high similarities anyway if the sign of
     * enrichment is different.
     * 
     * -> This is done to prevent early elimination of terms which can't be
     *    considered redundant as they have different contexts.
     * -> It also eliminates the awkward situation of having to compare
     *    magnitudes of positive and negative properties (i.e. is a 2-fold
     *    enrichment less important than a 3-fold depletion??)
     */
    for ( Integer curKeyX : similMat.keySet() ) {
      
      HashMap<Integer,Double> curHash = similMat.get(curKeyX);
      
      for ( Integer curKeyY : curHash.keySet() ) {
        if (curKeyX.equals(curKeyY) ) continue;
        Double propX =
                termProperties.get( myGo.get(curKeyX), criterion );
        Double propY =
                termProperties.get( myGo.get(curKeyY), criterion );
        if ( propX == null || propY == null || Double.isNaN(propX) || Double.isNaN(propY) )
          continue;
        /*propX = Math.log( Math.max(propX, 1e-2) );
        propY = Math.log( Math.max(propY, 1e-2) );*/

        // GO categories with enrichment of different signs cannot eliminate one another
        if ( Math.signum(propX) != Math.signum(propY) ) 
          similMat.get(curKeyX).put(curKeyY, 0.0);
      }
      
    }  
    
    
    while ( true ) {
      // determine the maximum remaining pairwise similarity in the matrix
      maxRemainingSimil = -1.0; maxX = -1; maxY = -1;
      for ( Integer curKeyX : similMat.keySet() ) {
        HashMap<Integer,Double> curHash = similMat.get(curKeyX);
        for ( Integer curKeyY : curHash.keySet() ) {
          if (curKeyX.equals(curKeyY) ) continue;
          double curCorrel = curHash.get(curKeyY);
          if ( curCorrel > maxRemainingSimil ) {
            maxRemainingSimil = curCorrel;
            maxX = curKeyX; maxY = curKeyY;
          }
        }
      }
      if (maxRemainingSimil == -1.0) break; // this will end loop
      
      // now, either maxX or maxY have to be "removed" (i.e. marked as having
      // high dispensability, and all of their correlations in the matrix set 
      // to -1)


      int removeId;

      decideWhichToRemove: {

        // first condition - if one of the GO categories is "pinned" by the
        // user, it automatically wins the contest
        boolean pinnedX, pinnedY;
        pinnedX = termProperties.hasTermWithProperty(maxX, "pinned")
                && termProperties.get(maxX, "pinned") == 1.0;
        pinnedY = termProperties.hasTermWithProperty(maxY, "pinned")
                && termProperties.get(maxY, "pinned") == 1.0;
        if ( pinnedX && !pinnedY ) {
          removeId = maxY;
          break decideWhichToRemove;
        }
        if ( pinnedY && !pinnedX ) {
          removeId = maxX;
          break decideWhichToRemove;
        }

        // second condition - if one GO category is 'very general' and the other
        // GO category is not, the 'very general' category automatically loses
        if ( goTermSizes != null ) {
          double proportionX = goTermSizes.getFreqGuessIfUnknown(maxX,myGo);
          double proportionY = goTermSizes.getFreqGuessIfUnknown(maxY,myGo);
          if ( !Double.isNaN(proportionX) && !Double.isNaN(proportionY) ) {
            if ( proportionX > genegroupGeneralityThreshold &&
                    proportionY < genegroupGeneralityThreshold ) {
              removeId = maxX;
              break decideWhichToRemove;
            }
            else if ( proportionY > genegroupGeneralityThreshold &&
                    proportionX < genegroupGeneralityThreshold ) {
              removeId = maxY;
              break decideWhichToRemove;
            }
          }
        }

        // now, check the properties to see which term we like best
        Double propX = termProperties.get( myGo.get(maxX), criterion );
        Double propY = termProperties.get( myGo.get(maxY), criterion );

        if ( altCriterion != null &&
                (propX == null || Double.isNaN(propX) || propY == null || Double.isNaN(propY) ) ) {
          propX = termProperties.get( myGo.get(maxX), altCriterion );
          propY = termProperties.get( myGo.get(maxY), altCriterion );
        }

        if (propX == null || propY == null) {
          propX = 0.0; propY = 0.0;
        }

        /* AbsLogProp transformation is necessary for eliminating by enrichment;
         * this makes 0.001 more desirable than 0.01, while 100 is more desirable
         * by 10 (i.e. keeps the number further away from 1).
         *
         * AbsLogProp also doesn't hurt when eliminating by p-value, which is by
         * definition between 0 and 1, and the greater absolute logarithm is
         * always desirable: e.g. 1e-20 is better than 1e-10.
         */
        if ( absLogProp ) {
          propX = Math.abs( Math.log10( Math.max(propX, 1e-300) ) );
          propY = Math.abs( Math.log10( Math.max(propY, 1e-300) ) );
        }

        // if propX and propY are within the confidence interval, force equality
        boolean xWasGreaterBeforeForcedEquality = false;
        if ( Math.abs(propX-propY)/((propX+propY)/2) < confidenceInterval ) {
          xWasGreaterBeforeForcedEquality = ( propX > propY );
          propY = propX;
        }

        if ( propX > propY ) {
          removeId = keepGreater ? maxY : maxX;
        } else if ( propX < propY ) {
          removeId = keepGreater ? maxX : maxY;
        } else { // they are equal.. here the fun begins :)

          // check if the two are linked by a parent-child relationship
          Collection<GOTerm> parentsOfX = myGo.get(maxX).getAllParents();
          Collection<GOTerm> parentsOfY = myGo.get(maxY).getAllParents();
          // also get sizes of X and Y, we'll need this to decide which to keep
          double sizeOfX = 1, sizeOfY = 1;
          if ( goTermSizes != null ) {
            sizeOfX = goTermSizes.getSizeGuessIfUnknown( maxX, myGo );
            sizeOfY = goTermSizes.getSizeGuessIfUnknown( maxY, myGo );
          }

          if ( parentsOfX.contains( myGo.get(maxY) ) ) {

            /* Y is a parent of X. Now check if Y is constituted mostly of X
             * by comparing their log-sizes - if yes, keep X (the more
             * specific term). Otherwise, keep Y, the parent term.
             *
             * TODO: This should be reworked into something like "if X has most of
             * the flagged/non-flagged genes of Y, keep X". Flagged genes
             * would be more interesting in enriched categories, and non-flagged
             * in depleted ones. See example:
             *
             * GO:0015986	1674 flagged, 1865 not; p=-300; enr=3.795; ATP_synthesis_coupled_proton_transport
             * GO:0006818	1804 flagged, 3170 not; p=-300; enr=2.909; hydrogen_transport (parent)
             *
             * Currently, we keep the parent here, and we obviously shoud not.
             */
            if ( (sizeOfY - sizeOfX) / sizeOfY < 0.25 )
              removeId = maxY; // Y constituted mostly of X
            else
              removeId = maxX; // Y constituted only partly of X

          } else if ( parentsOfY.contains( myGo.get(maxX) ) ) {

            /* X is a parent of Y. Rule as above, only reversed. */
            if ( (sizeOfX - sizeOfY) / sizeOfX < 0.25 )
              removeId = maxX; // X constituted mostly of Y
            else
              removeId = maxY; // X constituted only partly of Y

          } else { // X and Y are not parent-child related - keep the more
                   // desirable value regardless of the confidence interval

            if (xWasGreaterBeforeForcedEquality)
              removeId = keepGreater ? maxY : maxX;
            else
              removeId = keepGreater ? maxX : maxY;

          }

        }

      } // end block "decideWhichToRemove:"


      int dispensedById = ( removeId == maxX ? maxY : maxX );

      termProperties.put(myGo.get(removeId), "dispensability", maxRemainingSimil);
      termProperties.put(myGo.get(removeId), "dispensedBy",  (double) dispensedById);

      for ( Integer curKeyX : similMat.keySet() ) {
        Map<Integer,Double> curHash = similMat.get(curKeyX);
        for ( Integer curKeyY : curHash.keySet() ) {
          if ( curKeyX.intValue()==removeId || curKeyY.intValue()==removeId )
            curHash.put(curKeyY, -1.0);
        }
      }
      
    } // end while ( maxRemainingCorrel > -1.0 )
            
            
  }
  
  

  /**
   * Makes a Weka Instances object by performing these steps:
   * <ul>
   *   <li> makes similarity matrix from list of GO terms and converts to Instances
   *   <li> computes 'dispensability' of each GO term (stored into GOTerm properties)
   *   <li> removes all dispensable GO Terms from similarity matrix (you can skip this step
   * by passing 1.0 as dispensabilityCutoff)
   *   <li> performs PCA on the matrix
   *   <li> adds GO ids and names to Instances
   *   <li> adds all stored properties to Instances
   * </ul>
   *
   * @param dispensabilityCriterion A name of the Genes' property used for
   * selection. Here, larger is assumed to be better; transform this value
   * appropriately before passing it to this function.
   *
   * @param dispensabilityCutoff Must be between 0.0 and 1.0 (although either
   * of these extreme values are not very sensible).
   *
   * @param measure Specify a measure to use (SimRel is commonly used).
   *
   * @return Weka Instances object.
   */
  public Instances makeOntologram( GoTermProperties termProperties, String dispensabilityCriterion,
          double dispensabilityCutoff, SemanticSimilarityScore measure, GoTermSizes goTermSizes ) throws Exception {

    if ( dispensabilityCutoff < 0.0 || dispensabilityCutoff > 1.0 )
      throw new IllegalArgumentException("Illegal value of " +
              "dispensabilityCutoff passed to makeOntologram()");

    this.trimHugeTermLists(termProperties, dispensabilityCriterion,
            goTermSizes);

    TermSimilarityMatrix myMatrix;
    /*if ( this.serverUrl != null ) {
      myMatrix = new TermSimilarityMatrix(termList,
              measure, myGo, this.serverUrl);
    } else if ( this.genePlex != null ) {
      myMatrix = new TermSimilarityMatrix( termList,
              genePlex, myGo);
    } else { */
      myMatrix = new TermSimilarityMatrix( termList,
              measure, myGo, this.termSizes);
    //}

    termList = (GoTermList) myMatrix.supportedTerms.clone();

    myMatrix.calculateUniqueness(termProperties);

    //System.err.println( myInst.toString() );
    Long millis;
    millis = System.currentTimeMillis();

    // this step adds/alters properties of GO Terms in the TermSimilarityMatrix
    // properties changed: "dispensability", "dispensedBy"
    // also, the "supportedTerms" list changes
    markDispensableTerms( termProperties, dispensabilityCriterion, 
            "uniqueness",
            myMatrix, goTermSizes );
    //System.out.println("Done with markDispensableTerms. Took "
    //       + (System.currentTimeMillis() - millis) + " millis.");

    millis = System.currentTimeMillis();
    // warning - after this, myMatrix.supportedTerms will only contain a subset
    // of terms in the original termList
    myMatrix.filterByProperty(termProperties, "dispensability", -Double.MAX_VALUE, dispensabilityCutoff);
    // these Instances have no class attribute and generally no other attributes
    // other than the pairwise distances themselves
    Instances myInst = myMatrix.toWekaInstances(description, 1.0);
    //Instances myInstRootDist = myMatrix.toWekaInstances(description, 0.25);
    //System.out.println("Done with filtering matrix and making Weka Instances. Took "
    //        + (System.currentTimeMillis() - millis) + " millis.");

    millis = System.currentTimeMillis();

    final int numPcs = 2;

    // skip PCA if we have 2 or 1 attributes/instances
    if ( myInst.numInstances() > numPcs  && myInst.numAttributes() > numPcs ) {
      try {
        //myInst = reduceDimensionality(myInst, numPcs);
        myInst = reduceDimensionalityMds(myInst, numPcs);
        //myInstRootDist = reduceDimensionality(myInstRootDist, 4);
      } catch  (WekaException ex) {
        /* Weka's PrincipalComponents used in reduceDimensionality() will 
         * delete attributes that have only one distinct value without asking,
         * prior to applying PCA.
         * This may (very rarely) cause a situation where all attributes are 
         * deleted, and a WekaException will be thrown. In these cases,
         * we'll proceed without PCA and use the original attributes for 
         * visualization. */
        System.err.printf( "While visualizing task: %s, reduceDimensionality()" +
                " has caused a WekaException:" +
                " %s. The original distance matrix will be used for output.",
                description, ex.toString());
      }
    }
    // let's scatter those eeevil overlapping points!
    //myInst = scatterPoints(myInst, 0, 1);
    
    /* For the first four PCs, add new properties to the GO Terms with their
     * values. */
    for ( int i = 0 ; i < myInst.numInstances() ; i++ ) {
      GOTerm curTerm = myMatrix.supportedTerms.get(i);  // order of terms in TermList should correspond to order of instances
      for ( int j = 1; j <= numPcs; j++ ) {
        if ( myInst.numAttributes()>=j ) { // PC is available (generally, it should be, except for very small distance matrices)
          termProperties.put(curTerm, "PC_"+j, myInst.instance(i).value(j-1));
        } else { // PC is not available (the distance matrix was very small)
          termProperties.put(curTerm, "PC_"+j, Math.random());
        }
      }
    }

    //System.out.println("Done with PCA. Took " + (System.currentTimeMillis() - millis) + " millis.");


    // In the weka.core.Instances object,
    // append one attribute for each Property present in all relevant GOTerm-s
    Set<String> allTermProperties = new LinkedHashSet<String> ();
    for ( GOTerm curGOTerm : myMatrix.supportedTerms )
      allTermProperties.addAll( termProperties.getAllProps(curGOTerm).keySet() );

    for ( String curProp : allTermProperties ){

      if ( curProp.equals(isTestBString) )
        continue;
      if ( curProp.contains("dist_0") )
        continue;
      if ( curProp.startsWith("PC_") )
        continue;

      Attribute newAtt = new Attribute( curProp );
      myInst.insertAttributeAt( newAtt, myInst.numAttributes() );
        
      for ( int i = 0 ; i < myInst.numInstances() ; i++ ) {

        Double val = termProperties.get(myMatrix.supportedTerms.get(i), curProp);
        if ( val == null ) { // this specific property is not defined for this specific GO Term
          myInst.instance(i).setMissing( myInst.numAttributes()-1 );
        } else { // the property is defined here
          if ( curProp.equals( pValueString ) )
            val = Math.log10( Math.max(val, 1e-300) );
          myInst.instance(i).setValue( myInst.numAttributes()-1, val);
        }

      }
      
    }


    // In the weka.core.Instances object, append before all other attributes
    // ... name of organism group (same as filename)
    // ... number of organisms in group, and total number of organisms where propery is defined, in form (xx|yy)
    // ... and the GO ID (string attribute)
    // and after all other attributs
    // ... the term name (also a string attribute) after all attributes

    myInst.insertAttributeAt(new Attribute("GeneGroup", (FastVector) null), 0);
    myInst.insertAttributeAt(new Attribute("OrgsIn|Total", (FastVector) null), 0);
    myInst.insertAttributeAt(new Attribute("%OrgGroup", (FastVector) null), 0);

    myInst.insertAttributeAt(new Attribute("GO_name", (FastVector) null), myInst.numAttributes() );

    for ( int i = 0 ; i < myInst.numInstances() ; i++ ) {
      
      myInst.instance(i).setValue( 0, description );
      myInst.instance(i).setValue( 1, orgsInTotal );
      myInst.instance(i).setValue( 2,
              String.format("GO:%07d", termList.get(i).getId()) );

      myInst.instance(i).setValue( myInst.numAttributes()-1, "GO:" + termList.get(i).getId() + "_" + termList.get(i).getSafeName() );
    }

    return myInst;
  }

  /**
   * Takes a list of GO terms from the current TermListVisualizer object,
   * checks if the list is 'huge' (i.e. above a certain threshold) and if it 
   * is, makes sure that it fits within that size by throwing out GO terms that:
   * 
   * (a) have 'bad' p-values/enrichments/whatever - note: the GoTermProperties
   *     must have been transformed prior to invoking TermListVisualizer so that
   *     larger is always better
   * 
   * OR, if the values above were not supplied by the user
   * 
   * (b) are very general.
   */
  private void trimHugeTermLists( GoTermProperties termProperties,
          String dispensabilityCriterion,
          GoTermSizes goTermSizes ) {

    if ( this.termList.size() <= maxAllowedGoListSize )
      return;

    boolean allTermsHaveProperty = true;
    ArrayList<Double> propValues = new ArrayList<Double>();
    for ( GOTerm curTerm : this.termList ) {
      if ( termProperties.hasTermWithProperty(curTerm, dispensabilityCriterion) && ! Double.isNaN(termProperties.get(curTerm, dispensabilityCriterion)) ) {
        propValues.add( termProperties.get(curTerm, dispensabilityCriterion));
      } else {
        allTermsHaveProperty = false; //uh-oh, this won't work; use generality instead
        break;
      }
    }

    if (!allTermsHaveProperty) {
      // resort to using Go Term generality instead
      propValues.clear();
      for ( GOTerm curTerm : this.termList ) {
        propValues.add( - goTermSizes.getFreqGuessIfUnknown(curTerm.getId(), curTerm.getGeneOntology()) );
      }
    }


    //Double[] arr = (Double[]) propValues.toArray();
    Collections.sort(propValues);
    //Arrays.sort(arr);
    double threshold = propValues.get(this.termList.size() - maxAllowedGoListSize ); // keep maxAllowedGoListSize biggest values

    Iterator<GOTerm> it = this.termList.iterator();
    while( it.hasNext() ){
      double valToTest;
      GOTerm curTerm = it.next();
      if ( allTermsHaveProperty ) {
        valToTest = termProperties.get(curTerm, dispensabilityCriterion);  // remove 'undesirable' GO categories, e.g. with large p-values
      } else {
        valToTest = - goTermSizes.getFreqGuessIfUnknown(curTerm.getId(), curTerm.getGeneOntology()); // remove general (i.e. not highly informative) GO categories
      }
      if (valToTest < threshold) {
        it.remove();
      }
    }

  }


  /** Lists larger than this are trimmed prior to visualization. */
  public static int maxAllowedGoListSize = 350;
  

  /**
   * A helper class containing two distinct StringBuilders (edges, nodes), and
   * a method to write their contents out to two files.
   */
  protected class CytoscapeInputFiles {
    
    StringBuilder edges = null;
    StringBuilder nodes = null;
    
    public void writeOutToFiles ( String edgesFilename, String nodesFilename )
            throws IOException {

      if ( edges == null || edges.length() == 0 )  // don't write empty files
        return;

      PrintWriter pw = new PrintWriter(new FileWriter(nodesFilename, false), true);
      pw.print( nodes );
      pw.close();
      pw = new PrintWriter(new FileWriter(edgesFilename, false), true);
      pw.print( edges );
      pw.close();
      
    }
    
  }
  

  public class Edge{
      int sourceID;
      int destinationID;
      HashMap<Object,Object> properties;
  }
  public class Node{
      int ID;
      HashMap<Object,Object> properties;
  }
  
  public class OntoloGraph{

    ArrayList<Edge> edges;
    ArrayList<Node> nodes;

    public void writeAdditionalAttributes(){
        double min = Double.POSITIVE_INFINITY;
        for (int i=0;i<nodes.size();++i)
            if ((Double)nodes.get(i).properties.get("value")<min)
                min=(Double)nodes.get(i).properties.get("value");
        double max = Double.NEGATIVE_INFINITY;
        for (int i=0;i<nodes.size();++i)
            if ((Double)nodes.get(i).properties.get("value")>max)
                max = (Double)nodes.get(i).properties.get("value");

        for (int i=0;i<nodes.size();++i){
            double value = (Double)nodes.get(i).properties.get("value");
            String color = "";
            if (value <0){
                color = Integer.toHexString(255-(int)(255*value/(min-0)));
                if (color.length()==1)
                    color = "0"+color;

                color = "#ff"+color+color;
            }
            else {
                color = Integer.toHexString(255 - (int)(255*value/(max-0)));
                if (color.length()==1)
                    color = "0"+color;
                color = "#"+color+"ff"+color;
            }
            nodes.get(i).properties.put("color", color);
        }


        /*min = Double.POSITIVE_INFINITY;
        for (int i=0;i<nodes.size();++i)
            if ((Double)nodes.get(i).properties.get("log_size")<min)
                min=(Double)nodes.get(i).properties.get("log_size");
        max = Double.NEGATIVE_INFINITY;
        for (int i=0;i<nodes.size();++i)
            if ((Double)nodes.get(i).properties.get("log_size")>max)
                max = (Double)nodes.get(i).properties.get("log_size");
        for (int i=0;i<nodes.size();++i){
            double value = (Double)nodes.get(i).properties.get("log_size");
            value = 20+130*value/(max-min);
            nodes.get(i).properties.put("size", value);
        }*/

        min = Double.POSITIVE_INFINITY;
        for (int i=0;i<edges.size();++i)
            if ((Double)edges.get(i).properties.get("similarity")<min)
                min=(Double)edges.get(i).properties.get("similarity");
        max = Double.NEGATIVE_INFINITY;
        for (int i=0;i<edges.size();++i)
            if ((Double)edges.get(i).properties.get("similarity")>max)
                max = (Double)edges.get(i).properties.get("similarity");
        for (int i=0;i<edges.size();++i){
            double value = (Double)edges.get(i).properties.get("similarity");
            value = 1+4*value/(max-min);
            edges.get(i).properties.put("thickness", value);
        }
    }

    /** Format used as input for Cytoscape Web. */
    public String writeOutToGraphML() throws IOException{
        StringBuilder fileContent = new StringBuilder(10000);
        Formatter formatter = new Formatter(fileContent, Locale.US);
        formatter.format("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
        formatter.format("<graphml xmlns=\"http://graphml.graphdrawing.org/xmlns\"\n");
        formatter.format("xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"\n");
        formatter.format("xsi:schemaLocation=\"http://graphml.graphdrawing.org/xmlns\n");
        formatter.format("http://graphml.graphdrawing.org/xmlns/1.0/graphml.xsd\">\n");
        formatter.format("<key id=\"description\" for=\"node\" attr.name=\"description\" attr.type=\"string\"/>\n");
        formatter.format("<key id=\"value\" for=\"node\" attr.name=\"value\" attr.type=\"float\"/>\n");
        formatter.format("<key id=\"color\" for=\"node\" attr.name=\"color\" attr.type=\"string\"/>\n");
        formatter.format("<key id=\"log_size\" for=\"node\" attr.name=\"log_size\" attr.type=\"float\"/>\n");
        formatter.format("<key id=\"weight\" for=\"edge\" attr.name=\"weight\" attr.type=\"float\"/>\n");
        formatter.format("<graph id=\"REViGOgraph\" edgedefault=\"undirected\">\n");

        for (int i=0;i<nodes.size();++i){
            formatter.format("<node id=\"GO:%d\">\n",nodes.get(i).ID);
            formatter.format("\t<data key=\"description\">%s</data>\n", nodes.get(i).properties.get("description").toString().replaceAll("'", "´"));
            formatter.format("\t<data key=\"value\">%4.3f</data>\n", nodes.get(i).properties.get("value"));
            formatter.format("\t<data key=\"color\">%s</data>\n", nodes.get(i).properties.get("color"));
            formatter.format("\t<data key=\"log_size\">%4.3f</data>\n", nodes.get(i).properties.get("log_size"));
            formatter.format("</node>\n");
        }
        for (int i=0;i<edges.size();++i) {
            formatter.format("<edge source=\"GO:%d\" target=\"GO:%d\">\n", edges.get(i).sourceID, edges.get(i).destinationID);
            formatter.format("\t<data key=\"weight\">%4.3f</data>\n", edges.get(i).properties.get("similarity"));
            formatter.format("</edge>\n");
        }
        formatter.format("</graph>\n</graphml>");

        return fileContent.toString();
    }

    /** Format used as input for Cytoscape (standalone or via Java Web start) */
    public String writeOutToXGMML() throws IOException{
        StringBuilder fileContent = new StringBuilder(10000);
        Formatter formatter = new Formatter(fileContent, Locale.US);

        formatter.format("<graph label=\"Cytoscape Web\" directed=\"0\" Graphic=\"1\" xmlns:dc=\"http://purl.org/dc/elements/1.1/\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" xmlns:rdf=\"http://www.w3.org/1999/02/22-rdf-syntax-ns#\" xmlns:cy=\"http://www.cytoscape.org\" xmlns=\"http://www.cs.rpi.edu/XGMML\">\n");
        formatter.format("<att name=\"documentVersion\" value=\"0.1\"/>\n");

        for (int i=0;i<nodes.size();++i){
            formatter.format("<node id=\"GO:%d\" label=\"%s\">\n",
                    nodes.get(i).ID,nodes.get(i).properties.get("description"));
            formatter.format("\t<att type=\"string\" name=\"description\" value=\"%s\"/>\n",
                    nodes.get(i).properties.get("description"));
            formatter.format("\t<att type=\"real\" name=\"value\" value=\"%4.3f\"/>\n",
                    nodes.get(i).properties.get("value"));
            formatter.format("\t<att type=\"real\" name=\"log_size\" value=\"%4.3f\"/>\n",
                    nodes.get(i).properties.get("log_size"));
            formatter.format("\t<att type=\"real\" name=\"uniqueness\" value=\"%4.3f\"/>\n",
                    nodes.get(i).properties.get("uniqueness"));
            formatter.format("\t<att type=\"real\" name=\"dispensability\" value=\"%4.3f\"/>\n",
                    nodes.get(i).properties.get("dispensability"));
            formatter.format("\t<graphics type=\"ELLIPSE\" x=\"%4.3f\" y=\"%4.3f\" fill=\"%s\"/>\n",
                    25*(Double)nodes.get(i).properties.get("PC_1"),
                    25*(Double)nodes.get(i).properties.get("PC_2"),
                    nodes.get(i).properties.get("color"));
            formatter.format("</node>\n");
        }
        for (int i=0;i<edges.size();++i){
            formatter.format("<edge id=\"e%d\" target=\"GO:%d\" source=\"GO:%d\" directed=\"false\" label=\"e104\">\n",
                    i,edges.get(i).sourceID,edges.get(i).destinationID);
            formatter.format("\t<att type=\"real\" name=\"weight\" value=\"%4.3f\"/>\n",
                    edges.get(i).properties.get("similarity"));
            formatter.format("<graphics width=\"%4.3f\" fill=\"#edecec\"/>",
                    (Double)edges.get(i).properties.get("thickness"));
            formatter.format("</edge>\n");
        }
        formatter.format("</graph>");

        return fileContent.toString();
    }
      
  }




  public OntoloGraph makeSimpleThresholdOntolograph(
        GoTermProperties termProps,
        Double dispensabilityCutoff,
        Double similarityThreshold,
        String dispensabilityCriterion,
        SemanticSimilarityScore measure,
        GoTermSizes goTermSizes)
        throws Exception {

    /* if ( !measure.equals(SemanticSimilarityScore.SIMREL) )
    throw new IllegalArgumentException("Currently only SIMREL is supported" +
            " as a semantic similarity measure in makeOntolograph (others" +
            " don't have null distributions determined)."); */
    OntoloGraph result = new OntoloGraph();
    result.edges = new ArrayList<Edge>();
    result.nodes = new ArrayList<Node>();

    //CytoscapeInputFiles result = new CytoscapeInputFiles();
    if (termList.size() == 0)
      return null; // empty if no input data

    TermSimilarityMatrix myMatrix;
    /*if ( this.serverUrl != null ) {
      myMatrix = new TermSimilarityMatrix(termList,
              measure, myGo, this.serverUrl);
    } else if ( this.genePlex != null ) {
      myMatrix = new TermSimilarityMatrix( termList,
              genePlex, myGo);
    } else { */
      myMatrix = new TermSimilarityMatrix( termList,
              measure, myGo, this.termSizes);
    //}
    termList = myMatrix.supportedTerms;

    List<GOTerm> outerTermsProcessed = new GoTermList();
    Set<Integer> termIdsWithConnections = new HashSet<Integer>();

    int noOfEdges = myMatrix.size()*(myMatrix.size()-1)/2;
    double[] similarities = new double[noOfEdges];

    int counter = 0;
    for (GOTerm outerGoTerm : termList) {
      int outerGoId = outerGoTerm.getId();
      outerTermsProcessed.add(outerGoTerm);
      for (GOTerm innerGoTerm : termList) {
        int innerGoId = innerGoTerm.getId();
        // skip connections to self, and two-way connections
        if ( innerGoId == outerGoId ) continue;
        if ( outerTermsProcessed.contains(innerGoTerm) ) continue;

        if (termProps.get(outerGoTerm, "pinned") != 1.0){
            if ( termProps.getAllProps(outerGoTerm).containsKey("dispensability")
                       &&
                       termProps.get(outerGoTerm, "dispensability") > dispensabilityCutoff )
                  continue;
        }
        if (termProps.get(innerGoTerm, "pinned") != 1.0){
                    if ( termProps.getAllProps(innerGoTerm).containsKey("dispensability")
                       &&
                       termProps.get(innerGoTerm, "dispensability") > dispensabilityCutoff )
                  continue;
        }


        double simil = myMatrix.get( outerGoId ).get( innerGoId );
        similarities[counter++] = simil;
      }
    }
    Arrays.sort(similarities);
    double threshold;
    if ( noOfEdges == 0 ) {
      threshold = Double.NaN; // in this case, it won't matter what this value is anyhow..
    } else {
      threshold = similarities[ Math.max( (int)(noOfEdges*similarityThreshold) - 1, 0) ];
    }

    outerTermsProcessed = new GoTermList();
    termIdsWithConnections = new HashSet<Integer>();

    for (GOTerm outerGoTerm : termList) {
      int outerGoId = outerGoTerm.getId();

      outerTermsProcessed.add(outerGoTerm);

      for (GOTerm innerGoTerm : termList) {
        int innerGoId = innerGoTerm.getId();

        // skip connections to self, and two-way connections
        if ( innerGoId == outerGoId ) continue;
        if ( outerTermsProcessed.contains(innerGoTerm) ) continue;

        double simil = myMatrix.get( outerGoId ).get( innerGoId );


        if (termProps.get(outerGoTerm, "pinned") != 1.0){
            if ( termProps.getAllProps(outerGoTerm).containsKey("dispensability")
                       &&
                       termProps.get(outerGoTerm, "dispensability") > dispensabilityCutoff )
                  continue;
        }

        if (termProps.get(innerGoTerm, "pinned") != 1.0){
                    if ( termProps.getAllProps(innerGoTerm).containsKey("dispensability")
                       &&
                       termProps.get(innerGoTerm, "dispensability") > dispensabilityCutoff )
                  continue;
        }


        if ( simil >= threshold ) {
            Edge edge = new Edge();
            edge.properties = new HashMap<Object, Object>();
            edge.sourceID = outerGoId;
            edge.destinationID = innerGoId;
            edge.properties.put("similarity", simil);
            result.edges.add(edge);
            termIdsWithConnections.add( innerGoId );
            termIdsWithConnections.add( outerGoId );
        }
      }
    }

    termList.iterator();

    // the following part is commented out because before this routine,
    // it is expected that makeOntolograph() was called, where the uniqueness
    // and dispensability were computed

    /* myMatrix.calculateUniqueness(termProps);
    // WARNING - this step adds properties to GO Terms
    //markDispensableTerms( termProps, dispensabilityCriterion, "uniqueness", myMatrix, goTermSizes);

    int numPcs = 2;
    Instances myInst = myMatrix.toWekaInstances(description, 1.0);
    //myInst = reduceDimensionality(myInst, numPcs);
    myInst = reduceDimensionalityMds(myInst, numPcs);
    for ( int i = 0 ; i < myInst.numInstances() ; i++ ) {
      GOTerm curTerm = myMatrix.supportedTerms.get(i);  // order of terms in TermList should correspond to order of instances
      for ( int j = 1; j <= numPcs; j++ ) {
        if ( myInst.numAttributes()>=j ) { // PC is available (generally, it should be, except for very small distance matrices)
          termProps.put(curTerm, "PC_"+j, myInst.instance(i).value(j-1));
        } else { // PC is not available (the distance matrix was very small)
          termProps.put(curTerm, "PC_"+j, Math.random());
        }
      }
    } */

    /* Create a file describing the nodes. Start by enumerating all properties
     * of GO terms, then iterate term by term and write out all its properties
     * into appropriate columns. */
    Set<String> allTermProperties = new LinkedHashSet<String> ();
    for ( GOTerm curGoTerm : termList )
      allTermProperties.addAll( termProps.getAllProps(curGoTerm).keySet() );

    for ( GOTerm curGoTerm : termList ) {
        if (termProps.get(curGoTerm, "pinned") != 1.0){
            if ( termProps.getAllProps(curGoTerm).containsKey("dispensability")
                       &&
                       termProps.get(curGoTerm, "dispensability") > dispensabilityCutoff )
                  continue;
        }

        Node node = new Node();
        node.properties = new HashMap<Object, Object>();
        node.ID=curGoTerm.getId();
        node.properties.put("description", curGoTerm.getAbbreviatedName());

        String[] props = {"value", "log_size", "PC_1", "PC_2", "dispensability", "uniqueness"};
        for (int i=0;i<props.length;++i){
            String curProp = props[i];
            if ( termProps.getAllProps(curGoTerm).containsKey(curProp) ) {
                double val = termProps.get(curGoTerm, curProp);
                if ( curProp.equals( pValueString ) )
                    val = Math.log10( Math.max(val, 1e-300) );
                    node.properties.put(curProp, val);
            } else
                node.properties.put(curProp, 0);
        }

        result.nodes.add(node);
    }

    result.writeAdditionalAttributes();
    return result;
  }


















  
  /** 
   * Makes a CytoscapeInputFiles object, which contains two StringBuilders:
   * 
   * - one contains information about the edges of the graph, where edge length
   *   is semantic similarity of two GO terms - if it's above the similarity
   *   expected by chance (see SemanticSimilarityScore class)
   * 
   * - the other describes the nodes in the graph (GO categories), where
   *   'dispensability' of each GO term is computed and outputted to this file,
   *   in addition to its other properties - enrichment, size, p-value...
   *
   * @param dispensabilityCriterion A name of the Genes' property used for
   * selection.
   * 
   * @param measure Specify a measure to use (SimRel is commonly used).*
   * @throws java.lang.Exception
   */
  protected CytoscapeInputFiles makeOntolograph(
          GoTermProperties termProps,
          String dispensabilityCriterion,
          SemanticSimilarityScore measure,
          GoTermSizes goTermSizes)
          throws Exception {
    
    if ( !measure.equals(SemanticSimilarityScore.SIMREL) )
      throw new IllegalArgumentException("Currently only SIMREL is supported" +
              " as a semantic similarity measure in makeOntolograph (others" +
              " don't have null distributions determined).");

    CytoscapeInputFiles result = new CytoscapeInputFiles();
    if (termList.size() == 0) return result; // empty if no input data

    TermSimilarityMatrix myMatrix;
    myMatrix = new TermSimilarityMatrix( termList,
            measure, myGo, this.termSizes);
    termList = myMatrix.supportedTerms;
    
    /* create a file for the graph structure - for each term-term
     * similarity which is above random, add an edge */
    result.edges = new StringBuilder(10000);
    Formatter fmEdges = new Formatter(result.edges, Locale.US);
    fmEdges.format("term1\tterm2\tsimilarity\n");
    
    List<GOTerm> outerTermsProcessed = new GoTermList();
    Set<Integer> termIdsWithConnections = new HashSet<Integer>();
    
    for (GOTerm outerGoTerm : termList) {
      int outerGoId = outerGoTerm.getId();
      
      outerTermsProcessed.add(outerGoTerm);
      
      for (GOTerm innerGoTerm : termList) {
        int innerGoId = innerGoTerm.getId();
        
        // skip connections to self, and two-way connections
        if ( innerGoId == outerGoId ) continue; 
        if ( outerTermsProcessed.contains(innerGoTerm) ) continue;

        double simil = myMatrix.get( outerGoId ).get( innerGoId );
        double threshold = measure.nullDistForNamespace(namespace, true)[80];
        
        if ( simil > threshold ) {
          fmEdges.format("GO:%d\tGO:%d\t%4.3f\n", outerGoId, innerGoId, simil);
          termIdsWithConnections.add( innerGoId );
          termIdsWithConnections.add( outerGoId );
        }
        
      }
    }

    /* There is a possibility that a term has no connections with the strength
     * above the threshold value; this means it would not appear in the
     * ".graphEdges" file, and that it would not be loaded into Cytoscape
     * at all. This is solved by giving those terms a single self-connection. */
    for (GOTerm curGoTerm : termList) {
      int curId = curGoTerm.getId();
      if ( ! termIdsWithConnections.contains(curId) )
        fmEdges.format("GO:%d\tGO:%d\t%4.3f\n", curId, curId, 1.0);
    }


    termList.iterator();

    myMatrix.calculateUniqueness(termProps);
    /* WARNING - this step adds properties to GO Terms */
    markDispensableTerms( termProps, dispensabilityCriterion, "uniqueness",
            myMatrix, goTermSizes);
    
    /* Create a file describing the nodes. Start by enumerating all properties
     * of GO terms, then iterate term by term and write out all its properties
     * into appropriate columns. */
    Set<String> allTermProperties = new LinkedHashSet<String> ();
    for ( GOTerm curGoTerm : termList )
      allTermProperties.addAll( termProps.getAllProps(curGoTerm).keySet() );

    result.nodes = new StringBuilder(10000);
    Formatter fmNodes = new Formatter(result.nodes, Locale.US);

    fmNodes.format("term\tdescription"); // first write header
    for ( String curProp : allTermProperties )
      fmNodes.format("\t%s", curProp);
    fmNodes.format("\n");
    
    for ( GOTerm curGoTerm : termList ) {
      fmNodes.format("GO:%d\t%s", curGoTerm.getId(), curGoTerm.getAbbreviatedName() );
      for ( String curProp : allTermProperties ) {
        if ( termProps.getAllProps(curGoTerm).containsKey(curProp) ) {
          double val = termProps.get(curGoTerm, curProp);
          if ( curProp.equals( pValueString ) )
            val = Math.log10( Math.max(val, 1e-300) );
          fmNodes.format( "\t%4.3f", val  );
        } else
          fmNodes.format("\t?");
      }
      fmNodes.format("\n");
    }
    
    return result;
    
  }
  
  
  
  /**
   * Demo of the TermListVisualizer.
   */
  private static void demo( ) throws Exception {
    
    GeneOntology myGO = new GeneOntology();
    
    GoTermList myTermList = new GoTermList ();
    myTermList.add(myGO.get(27));
    GoTermProperties termProps = new GoTermProperties(myGO);
    termProps.put(myGO.get(27), enrichString, 3.563);
    termProps.put(myGO.get(27), log10SizeString, 100.0);
    termProps.put(myGO.get(27), pValueString, 4.5e-20);
    myTermList.add(myGO.get(16071));
    termProps.put(myGO.get(16071), enrichString, 3.319);
    termProps.put(myGO.get(16071), log10SizeString, 52.0);
    termProps.put(myGO.get(16071), pValueString, 2.5e-10);
    myTermList.add(myGO.get(6183));
    termProps.put(myGO.get(6183), enrichString, 3.2);
    termProps.put(myGO.get(6183), log10SizeString, 101.0);
    termProps.put(myGO.get(6183), pValueString, 1.5e-10);
    myTermList.add(myGO.get(6228));
    termProps.put(myGO.get(6228), enrichString, 3.2);
    termProps.put(myGO.get(6228), log10SizeString, 101.0);
    termProps.put(myGO.get(6228), pValueString, 1.8e-30);
    myTermList.add(myGO.get(6241));
    termProps.put(myGO.get(6241), enrichString, 3.2);
    termProps.put(myGO.get(6241), log10SizeString, 101.0);
    termProps.put(myGO.get(6241), pValueString, 1.6e-10);
    myTermList.add(myGO.get(51262));
    termProps.put(myGO.get(51262), enrichString, 3.033);
    termProps.put(myGO.get(51262), log10SizeString, 46.0);
    termProps.put(myGO.get(51262), pValueString, 1.5e-20);
    myTermList.add(myGO.get(6796));
    termProps.put(myGO.get(6796), enrichString, 2.971);
    termProps.put(myGO.get(6796), log10SizeString, 89.0);
    termProps.put(myGO.get(6796), pValueString, 1.5e-30);
    myTermList.add(myGO.get(6685));
    termProps.put(myGO.get(6685), enrichString, 2.965);
    termProps.put(myGO.get(6685), log10SizeString, 26.0);
    termProps.put(myGO.get(6685), pValueString, 1.5e-10);
    myTermList.add(myGO.get(6450));
    termProps.put(myGO.get(6450), enrichString, 2.953);
    termProps.put(myGO.get(6450), log10SizeString, 97.0);
    termProps.put(myGO.get(6450), pValueString, 1.5e-10);
    myTermList.add(myGO.get(6662));
    termProps.put(myGO.get(6662), enrichString, 2.873);
    termProps.put(myGO.get(6662), log10SizeString, 188.0);
    termProps.put(myGO.get(6662), pValueString, 1.5e-10);
    myTermList.add(myGO.get(55085));
    termProps.put(myGO.get(55085), enrichString, 0.286);
    termProps.put(myGO.get(55085), log10SizeString, 113.0);
    termProps.put(myGO.get(55085), pValueString, 1.5e-10);
    myTermList.add(myGO.get(34038));
    termProps.put(myGO.get(34038), enrichString, 0.2593);
    termProps.put(myGO.get(34038), log10SizeString, 170.0);
    termProps.put(myGO.get(34038), pValueString, 1.5e-10);

    
    TermListVisualizer myGOTermListVis =
            new TermListVisualizer(myTermList, myGO,
            GoTermSizes.deserializeFromFile( new File(
            "Z:/databases/EBI_GOA_mappings/goa_2009-11-23/" +
            "goCatSizes_swissProtOnly/0.ser") ) );

    /*TermListVisualizer myGOTermListVis =
            new TermListVisualizer(myTermList, myGO, "http://funsimmat.bioinf.mpi-inf.mpg.de/xmlrpc.php" );*/


    Instances myInst = myGOTermListVis.makeOntologram(termProps,
            enrichString, 1.0, SemanticSimilarityScore.SIMREL, null);
    System.out.println(myInst);
           
    
  }



  /** are we looking for higher enrichments/depletions, or stronger p-values? */
  final static String propertyOfInterest = pValueString; //enrichString; // pValueString;



  /**
   * Takes the output file of GeneMonkey, with each row containing results
   * from a single enrichment test.
   * 
   * Produces many output files (one per OrganismGroup) with same results
   * 'located' in a space defined by a Gene Ontology term similarity metric,
   * and ranked by their 'dispensability' - see JavaDoc of
   * markDispensableTerms() for a definition.
   * 
   * @param args A single command line parameter is expected, giving the 
   * location of the GeneMonkey output file.
   * @throws java.lang.Exception
   */
  public static void main( String[] args ) throws Exception {
    
    if ( args.length == 0 ) {
      System.out.println("You didn't supply the name of the GeneMonkey output " +
              "file. The demo() function will be run instead.");
      demo();
      System.exit(0);
    }

    GeneOntology myGo = new GeneOntology();
    BufferedReader in = new BufferedReader(new FileReader(args[0]));
    String line;
    String curOrgGroup = null;
    GoTermSizes goSizes = GoTermSizes.deserializeFromResources( "goa_2009-11-23", 0 );

    /* For the current OrgGroup, stores as a string the number of organisms in
     * it and the total number of organisms where the Genome Propery (on
     * which this subdivision is based) is defined. Format: "11|135" */
    String orgsInTotal = null;

    //Connection conn = new hr.irb.mirrorBird.buildDatabase
    //        .GenomesDatabase().getConnection();
    //List<Integer> allGeneIds = hr.irb.mirrorBird.buildDatabase
    //        .GenomesDatabase.getSampleOfGeneIds(1000000);
    //List<Integer> allGeneIds = hr.irb.mirrorBird.buildDatabase
    //        .GenomesDatabase.getAllGeneIds();
    //GenePlex genePlex = new GenePlex(  // GenePlex used only for internal measure of similarity
    //        allGeneIds, true, myGo );  // no Genome Properties required


    /* the GeneGroups that are GO categories (having IDs<0) are separated by  
     the root node they belong to (molecular function, biological process and
     cellular component) and saved in three ArrayLists. Size and
     enrichment of GO categories are saved in the Properties HashMap for 
     each GO term. */
    GoTermList molFunList = new GoTermList();
    GoTermList biolProcList = new GoTermList();
    GoTermList celCompList = new GoTermList();
    GoTermList allNamespacesList = new GoTermList();
    GoTermProperties termProps = new GoTermProperties(myGo);
    
    while (true) { // ... test by test

      line = in.readLine();
      String[] parsedLine = null;

      if ( line != null ) {

        line = line.trim();

        if ( (line.equals("") || line.startsWith("%")) && in.ready() )
          continue;  // skip empty lines and comments / header

        parsedLine = line.split("\t");
        
      }
      
      // write information when we change OrganismGroups, or reach end of file
      // TODO fix bug - skipping of the last line
      if ( ( curOrgGroup != null && parsedLine != null && !parsedLine[0].equals(curOrgGroup) )
             || line == null ) {
        
        String filename = curOrgGroup; /* "allResults" ;*/

        CSVSaver saver = new CSVSaver();
         
        if (molFunList.size() > 0) {
          saver.setInstances(new TermListVisualizer( molFunList,
                  myGo,
                  goSizes,
                  /*"molecularFunction_" +*/ filename, orgsInTotal
                  ).makeOntologram(termProps, propertyOfInterest, 1.0, SemanticSimilarityScore.SIMREL, goSizes));
          saver.setFile( new File(filename + "_molecularFunction.csv") );
          saver.writeBatch();          
        }
        if (biolProcList.size() > 0) {
          saver.setInstances(new TermListVisualizer( biolProcList,
                  myGo,
                  goSizes,
                  /*"biologicalProcess_" +*/ filename, orgsInTotal
                  ).makeOntologram(termProps, propertyOfInterest, 1.0, SemanticSimilarityScore.SIMREL, goSizes));
          saver.setFile( new File(filename + "_biologicalProcess.csv") );
          saver.writeBatch();
        }
        if (celCompList.size() > 0) {
          saver.setInstances(new TermListVisualizer( celCompList,
                  myGo,
                  goSizes,
                  /*"cellularComponent_" +*/ filename, orgsInTotal
                  ).makeOntologram(termProps, propertyOfInterest, 1.0, SemanticSimilarityScore.SIMREL, goSizes));
          saver.setFile( new File(filename + "_cellularComponent.csv") );
          saver.writeBatch();

        }
        /*
         if (allNamespacesList.size() > 0) {
          saver.setInstances(new TermListVisualizer( allNamespacesList, 
                  myGo, goSizes, orgsInTotal
                  ).makeOntologram(termProps, propertyOfInterest, 1.0, SemanticSimilarityScore.SIMREL));
          saver.setFile( new File(filename + "_allNamespaces.csv") );
          saver.writeBatch();
        }
        */
        

        /*
        TermListVisualizer vis; 
        
        vis = new TermListVisualizer(biolProcList, myGo, "biologicalProcess_"+filename);
        vis.makeOntolograph(propertyOfInterest).writeOutToFiles(
                filename + "_biologicalProcess.graphEdges", 
                filename + "_biologicalProcess.graphNodes");

        vis = new TermListVisualizer(celCompList, myGo, "cellularComponent_"+filename);
        vis.makeOntolograph(propertyOfInterest).writeOutToFiles(
                filename + "_cellularComponent.graphEdges", 
                filename + "_cellularComponent.graphNodes");
        
        vis = new TermListVisualizer(molFunList, myGo, "molecularFunction_"+filename);
        vis.makeOntolograph(propertyOfInterest).writeOutToFiles(
                filename + "_molecularFunction.graphEdges", 
                filename + "_molecularFunction.graphNodes");

        vis = new TermListVisualizer(allNamespacesList, myGo, "allNamespaces_"+filename);
        vis.makeOntolograph(propertyOfInterest).writeOutToFiles(
                filename + "_allNamespaces.graphEdges", 
                filename + "_allNamespaces.graphNodes");
        */
        
        
        // next organism group!
        molFunList = new GoTermList();
        biolProcList = new GoTermList();
        celCompList = new GoTermList();
        allNamespacesList = new GoTermList();
        termProps = new GoTermProperties(myGo);
        //myGo.clearAllProperties();
        
      }

      if ( line == null )
        break; // end of file reached

      /* parsedLine[0] is the OrganismGroup - nicely formatted, without
       *   characters illegal in filenames
       * parsedLine[2] is the GeneGroup prefixed ID. It must be in form
       *   of "GO:0003775".
       * parsedLine[5] + parsedLine[6] is the GeneGroup size
       * parsedLine[9] is the statistical significance (p-value).
       * parsedLine[11] is the enrichment factor
       * parsedLine[12] is the description of GeneGroup.
       */

      curOrgGroup = parsedLine[0];
      orgsInTotal = parsedLine[1];

      /* skip the "custom gene groups" which can be created in GeneMonkey -
       * those are in fact not Gene Ontology categories and consequently have
       * no semantic similarity to other GO cats */
      if ( parsedLine[2].startsWith("GO:1") )
        continue;

      if ( parsedLine[2].startsWith("COG:") )  // we don't care about COGs
        continue;

      int curGoId = parsePrefixedId(parsedLine[2]);
      // int curGoId = parsePrefixedId(parsedLine[0]);

      if ( curGoId >= 0 ) // only Gene Ontology categories are accepted, and they
        continue;         // have negative IDs
                  
      curGoId = Math.abs(curGoId);
      
      if ( !myGo.containsKey(curGoId) ) { // a really weird thing, skip this
        System.err.println("GO ID not reconized by the GeneOntology object: "
                + curGoId);
        continue;
      }

      GOTerm curGOTerm = myGo.get(curGoId);

      termProps.put(curGOTerm, isTestBString,
              parsedLine[4].trim().equalsIgnoreCase("B") ? 1.0 : 0.0
              );

      termProps.put(curGOTerm, log10SizeString,
              Math.log10( Double.parseDouble(parsedLine[5].trim())
              + Double.parseDouble(parsedLine[6].trim()) )
              /* Math.log10( Double.parseDouble(parsedLine[2].trim()) ) */
              );
      
      termProps.put(curGOTerm, inGroupFlaggedCntString,
              Double.parseDouble(parsedLine[5].trim()) );
      termProps.put(curGOTerm, inGroupNoFlagCntString,
              Double.parseDouble(parsedLine[6].trim()) );
      termProps.put(curGOTerm, outOfGroupFlaggedCntString,
              Double.parseDouble(parsedLine[7].trim()) );
      termProps.put(curGOTerm, outOfGroupNoFlagCntString,
              Double.parseDouble(parsedLine[8].trim()) );
      termProps.put(curGOTerm, pValueString,
               Double.parseDouble(parsedLine[9].trim() ));
      termProps.put(curGOTerm, flagFreqString,
              Double.parseDouble(parsedLine[10].trim()) );

      termProps.put(curGOTerm, enrichString,
              Double.parseDouble(parsedLine[ 11 /*1*/ ].trim()) );

      if (curGOTerm.isBiologicalProcess()) {
        biolProcList.add(curGOTerm);
        termProps.put(curGOTerm, "namespace", 8150.0);
      } else if (curGOTerm.isCellularComponent()) {
        celCompList.add(curGOTerm);
        termProps.put(curGOTerm, "namespace", 5575.0);
      } else {
        molFunList.add(curGOTerm);
        termProps.put(curGOTerm, "namespace", 3674.0);
      }
      allNamespacesList.add(curGOTerm);
      
    }  
    
  }
    
 
  
  /**
   * Returns integer id from a string such as "GO:0003775" or "OMA:9871"
   * or "AllGenes". GO categories will have negative ids, OMA groups will have
   * positive IDs, and AllGenes has ID 0.
   */ 
  public static int parsePrefixedId( String prefixedId ) {
    if ( prefixedId.equalsIgnoreCase("AllGenes") )
      return 0;
    else if ( prefixedId.startsWith("GO:") )
      return -Integer.parseInt( prefixedId.substring(3) );
    else if ( prefixedId.startsWith("OMA:") )
      return -Integer.parseInt( prefixedId.substring(4) );
    else 
      throw new IllegalArgumentException("The prefixedId must either equal " +
              "\"AllGenes\", or start with \"GO:\" or \"OMA:\".");
  }  
  

  
  
  
  
  
  
}
